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We study the formation of large-scale coherent structures (a condensate) for a system of two weakly 
interacting classical waves. Using the coupled defocusing nonlinear Schrodinger (NLS) equations 
as a representative model, we focus on condensation in the phase mixing regime. We employ weak 
turbulence theory to provide a complete thermodynamic description of the classical condensation 
process. We show that the temperature and the condensate mass fractions are fully determined by 
the total number of particles in each component and the initial total energy. Moreover, we find 
that, at higher energies, condensation can occur in only one component. The theory presented 
provides excellent agreement with results of numerical simulations obtained by directly integrating 
the dynamical model. 
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Many classical systems in nature reveal the emergence 
of large scale coherent structures from a background 
irregular field characterised by small-scale fluctuations. 
Examples of systems that exhibit such behaviour include 
classical turbulence, nonlinear optics, superfluids, ultra- 
cold gases and Bose-Einstein Condensates (BECs), and 
the formation of the early universe. In certain regions 
of the parameter space, a large sub-class of these sys- 
tems can be described by a system of weakly nonlinear 
dispersive waves. A universal equation that governs the 
evolving field in such scenarios is then given by the Non- 
linear Schrodinger (NLS) equation. The process of self- 
organisation in the focusing NLS equation has been stud- 
ied in lj. It was found that a large-scale solitary wave 
tends to emerge from a sea of small-scale turbulent fluc- 
tuations. In this work, we concentrate on the defocusing 
NLS equation. This equation has been receiving increas- 
ing attention due to the experimental advances in BECs. 
In this context, the defocusing NLS equation corresponds 
to the Gross-Pitaevskii (GP) equation of a homoge- 
neous Bose gas. The GP equation has long been used 
as a model of a weakly interacting Bose gas at zero tem- 
perature. More recently, it has been argued Q that the 
GP equation can be used to model the long wavelength 
part of the spectrum of a BEC at finite temperatures. 
Numerical simulations conducted within this framework 
[H, [H have indeed confirmed this, revealing the ability 
of the model to capture the formation of a condensate 
from an initially turbulent state. 

With the rapid developments being made in experi- 
mental techniques, it is now possible to realize multi- 
component BECs formed by the simultaneous trapping 
and cooling of atoms in distinct spin or hypcrfine levels 
Q or different atomic species Q • The finite temperature 
dynamics of such Bose gas mictures is then goverened 
by a system of coupled NLS equations. While such sys- 
tems are of of interest in their own right, they also serve 
as idealised models to study symmetry-breaking phase 



transitions that are believed to have occurred in the early 
evolution of the universe. A specific example is given by 
the Kibble-Zurek mechanism [8[ of the formation of topo- 
logical defects following the rapid quench of the system 
below the point of second-order phase transitions. This 
scenario would correspond to the formation of cosmolog- 
ical vortons and springs that are analogous to the vortex 
ring-slaved wave and vortex ring-vortex ring complexes 
of BECs In addition to these physical examples, the 
coupled NLS equations are also encountered in the study 
of optical fibres and electromagnetic waves [icj |. Given 
the universality of these equations in the nonlinear sci- 
ences, an accurate thermodynamic description of the con- 
densation process in such a system could have significant 
implications in various branches of physics. 

In this Letter, we will generalise the results presented 
in ll|, that describe the condensation process of a one- 
component system, to a two-component system. We note 
that two-component systems tend to show a broad class 
of qualitatively different behavior depending on the rel- 
ative strengths of the intcrcomponcnt and intracompo- 
nent interactions. This can lead to contrasting regimes 
of condensation: the phase mixing regime and the phase 
separation regime [12j |. Consistent with our original as- 
sumptions of weakly interacting systems, we will focus 
exclusively on the phase mixing regime. 

We begin by considering the scenario of a system of 
two weakly interacting classical waves (e.g. Bose gases 
consisting of the same atoms in different spins) that have 
been rapidly cooled below the transition temperature. 
Their evolution from the resulting strongly nonequilib- 
rium initial state is then described by the coupled NLS 
equations given by 
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where ipi an( i ^2 are complex-valued classical fields cor- 
responding to each component, and a is the intracompo- 
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nent coupling constant. For the phase mixing regime, we 
require < a < 1. The dynamics governed by the above 
equations will conserve the total mass (number of parti- 
cles) given by N% = J\ipi\ 2 dx and N 2 = J \tp 2 \ 2 dx. In 
addition, the total energy (Hamiltonian) of the coupled 
system 
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will be conserved. Without loss of generality we shall 
assume that N% < N±. 

Despite the formal reversibility of the above Hamil- 
tonian system, the evolution of the nonlinear waves tpi 
and ip2 is nonintegrable giving rise to an effective dif- 
fusion in phase space. This results in an irreversible 
evolution to thermal equilibrium. By invoking the ran- 
dom phase approximation (assumption of quasi-Gaussian 
statistics), it is possible to derive closed irreversible ki- 
netic equations that describe the evolution of the sys- 
tem using Weak Turbulence Theory (WTT) 0. For a 
homogeneous system, we accomplish this by expressing 
the order parameters in terms of their Fourier transforms 
Vi = J^m /a k (i)e ik - x dk, = (2^372 J b k (t)e^dk. 
Substituting into Eq. we can derive expressions for 
the spectral number densities (a^Ok ) = n iS(ki — k?); 
(^ki^k 2 ) = — ka). Provided the nonlinearity in the 

system is sufficiently weak (i.e. N±/V <C 1; N 2 /V -c 1; 
a <C 1, where V is the volume of the system), we can 
derive the kinetic equation 
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xS(k + ki - k 2 - k 3 )<5(fc 2 + k\ - k\ - fc|)dkidk 2 dk 3 . 

Similarly, an equation for l k can be derived and follows 
directly from above by the symmetry in Eq. (JJJ). These 
equations conserve N\ = V J n k (t)dk, N 2 = V J l^(t)dk, 
and the kinetic energies Ei = V J fc 2 nk(i)<ik, E 2 = 
V J k 2 lk{t)d\t of each component. They admit two for- 
mal equilibrium solutions; the first corresponding to a 
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given by the Rayleigh- Jeans (RJ) distribution 
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Here, T is the thermodynamic temperature, and \i\ and 
\x 2 are the chemical potentials. Equation ([3]) satisfies a 
H-theorem for entropy growth which implies that the RJ 
distribution will be realized in practice. However, Eq. 
(f4|) provides only a formal solution since it leads to non- 
convergent expressions for Ni, N 2 , and the kinetic ener- 
gies Ei, E 2 as k — > 00. We recall that, for BECs, Eq. (TIJ 



is valid in the limit of large occupation numbers where 
a semi-classical description is valid. When ~ 1 and 
ik ~ lj Eq. (fT]) begins to break down and a full quantum 
mechanical treatment of the problem becomes necessary. 
To regularise the ultra-violet catastrophe, we introduce 
a cut-off k c such that n c i{\k c \) > 1, l cq {\k c \) > 1. This 
cut-off does not affect the equilibrium state provided a 
sufficiently large number of modes can be represented 
classically 0, [H[ ■ The reason is that a full quantum me- 
chanical description corresponds to a grand-canonical en- 
semble with fluctuations in particle number and energy. 
However, for sufficiently many modes, such fluctuations 
will be small and we can introduce the above truncation 
to reduce the system to a microcanonical ensemble where 
the number of particles and the energy are conserved. 

The RJ distributions corresponding to Eq. §3§ are only 
valid at sufficiently high energies when no condensate 
is present. At sufficiently low energies, Eq. ([3]) breaks 
down very rapidly giving way to the formation of a con- 
densate as elucidated in numerical simulations for a one- 
component system 0, HI EH and a two-component sys- 
tem (IBj . In the simplest scenario, condensates with zero 
wavenumbers are formed and are associated with the uni- 
form states provided we are in the phase mixing regime. 
If the condensates that form are strong in the sense that 
(Ni - n )/Ni < 1, {N 2 - l a )/N 2 « 1 (no = M 2 ; 
lo = \b \ 2 are the occupation numbers of the conden- 
sates in components 1 and 2, respectively), one can de- 
scribe the nonlinear dynamics at these later times by 
considering the evolution of small quasiparticle pertur- 
bations around the condensates. For our two-component 
system, we accomplish this by introducing the ansatz 
a k = [^r o 5{k)+a k {t)]e- m ^, b k = [VT 5(k)+b k (t)}e- u ^ . 
Upon substituting these expressions into the Fourier- 
transform representation of Eq. ([2|) , we introduce a trans- 
formation to diagonalize_terms in the Hamiltonian that 
are quadratic in a k and b k . A generalisation of Bogoli- 
ubov's transformation to diagonalise the Hamilto- 
nian in two-component systems was given in [l7| . To this 
end, we introduce the canonical variables Af = (£k, ?7k) T 
which are related to the original variables A = (ak, &k) T 
through the relation 
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where AP = (£ k >»?-k) r . and A = (a*_ k ,b*_ k ) 
preserve the properties of the Poisson bracket in the 
new basis, the transformation must satisfy the condition 
U + U +T — U~U~ T = I. A transformation that satisfies 
this condition and diagonaliscs the quadratic term is ob- 
tained when the elements uf, of the 2x2 transformation 
matrices U are given by 
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= \J k 2 JD^ denotes the ratios of the dis 
persion relations, where il = 
= [n + l a ± y/(n 
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ay/l /n . The resulting ex- 
pression for the Hamiltonian leads to kinetic equations 
for the canonical (quasiparticlc) densities (£ kl £ k ) = 
niS(ki — k 2 ); (vkiVk*,) = hS{ki - k 2 ) which are given by 
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where A£{ 2 = C^{ 2 5(k — ki — k 2 ) and C^> denote coef- 
ficients that will, in general, depend on ufj, n , l Q and 
a. Since we are only interested in equilibrium solutions, 
their precise form is not too important. The equation 
for ?k follows by symmetry of the dynamical equation. 
These kinetic equations are now given by resonant three- 
wave interactions and have a one parameter family of 
solutions given by n k q = T/Q+{k), and £ q = T/Q-(k), 
respectively. The condensates, therefore, strongly affect 
the equilibrium distributions of the quasiparticles. 

Using these equilibrium solutions, we can now derive a 
relation between the occupation numbers n Q , l a , and the 
total number of particles Ni, N 2 and the energy H. For 
a finite sized system, the Hamiltonian can be expressed 
in terms of the Fourier series %pi = J2 k a k exp(ik • x), 
ip 2 = -j= b k exp(ik • x). The Hamiltonian can then 
be written as H = H + H 2 + H 3 + H4 depending on 
how a a = a k =o and b = b k=0 , and non-zero modes en- 



ter the expansion: H a 



2V} 



2|a | 2 (7V 1 



|a | 2 )+2|6 | 2 (iV 2 -|fe | 2 )] + f [\a \ 2 N 2 + \b a \*Ni] , H 2 = 

Ekl(k 2 + + M a l< a -* + c - c -)] + ELP 2 + 

+ MttKb-K + cc.)] + EJfK6c<6*_ k + 

a 6> k 6 k + cc.)], H 3 = J2k 1 ,k 2 -k 3 iw( 2a ° ak i a k 2 a k 3 + 
2b b kl b* k2 b* k3 + cc.) + f{a b kl a* k2 b* k3 + b a kl a*k 2 K 3 + 

cc)](5 kl _ k2 _ k3 , H 4 = Efc 1 ,fc 2 ,fc 3 ,fc 4 [2T( afc i afc 2 a fc 3 a fe 4 + 



J2k 



de- 



hth^blj + f-(a fel 6fc 2 4 3 6^)]5 kl+k2 _k 3 -k4 
notes summation over k but excluding the k = mode. 
To relate the equilibrium distributions obtained from the 
kinetic equation (JT]) to this form of the Hamiltonian, we 



must diagonalise the quadratic term H 2 . Rewriting the 
Hamiltonian in terms of the basis N = (£ k ,77k) T , the 
quadratic part takes the form H 2 = Efc(^ + (^)£fc£fc + 
ft~ (k)rjkr) k ) . Now ensemble averaging the equations and 
using the equilibrium distributions n k q , and Z k q given 
above, we can express the occupation numbers of the 
two gases in the new basis as 
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The ensemble averaged Hamiltonian will have contribu- 
tions from only H , H 2 , and H4. Rewriting the re- 
sulting expression in the new basis, we obtain: (H) = 

Eo + E' fe (^i(fcR q + ^W) = Eo + 2TJ2' k 1. E - 
^ [N 2 + (JVi - n Q f + N 2 + (N 2 - Q 2 ] + f [NiN 2 ] denotes 
the energy of the ground state. Using either Eq. ([5]) or 
we can eliminate the temperature in the expression 
H) to obtain 
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This equation provides two algebraic relations for the two 
unknowns n a and l given H, Ni, N 2 . At intermedi- 
ate energies, we find that only one component condenses 
and the theory presented above breaks down (in prac- 
tice resulting in negative values of l ). In this part of 
the parameter space, we need to introduce the ansatz 
afc = [\/noS(k) + a k (t)]e~" lot only for the first compo- 
nent whilst assuming a purely continuous spectrum for 
the second component. This results in a coupled sys- 
tem of kinetic equations that are similar to Eqs. and 
([7]) above but where the first component is governed by 
three-wave resonances whereas the second is governed 
by four- wave resonances. These equations will have the 



equilibrium distributions n k 



and h 



where 



n k is now defined as above but with l D set to zero and 
cub = v + 2n k 2 /V is the classical single-component 
Bogoliubov dispersion relation. Therefore, when only 
one component is condensed, we obtain the expressions: 
Ni n = T£' fe and N 2 = The 

ensemble averaged Hamiltonian then takes the form 
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and now E = ±[N? + (JVi - n a ) 2 + 2N'$ + 2aN 1 N 2 }. 
This equation provides two algebraic expressions for the 
two unknowns n Q and fi2- We note that at the critical 
energy where the condensate in the second component 
vanishes, we have l = and \ii = 0. At this point, 
it can be shown that the equilibrium distributions given 
by 

n+(k) anc ^ n-(k) recm ce to and ^ and the two 

solutions given above match at the critical energy. This 
provides a solution for the thermodynamic state that is 
uniformly valid over the entire range of energies. 

To verify the theory, we numerically solved the coupled 
system in a cubic region with periodic boundary condi- 
tions using a pseudo-spectral method with spatial resolu- 
tion containing (64 3 ) grid points and a time-step of 0.01. 
The initial conditions were set using the random phase 
approximation (see e.g. [4]). To determine the properties 
of the system at equilibrium, we assumed that the er- 
godic hypothesis applies and used time-averages from our 
simulations to represent ensemble averages that arise in 
the theory presented above. Figure [1] presents results for 
the variation of the condensate mass fractions with (H) . 
The results shown reveal remarkable agreement with the 
theory presented over the entire range of energies. For 
very large energies where only a small condensate mass 
fraction of the first component is present, the theory de- 
viates very slightly from the results of the simulations. 
This occurs since many particles are non-condensed in 
this region violating the assumption of a strong conden- 
sate that is required for the theory. 

Given the excellent agreement between the theory and 
predictions, we can use Eq. (|lip to determine how the 
critical energy ((H) clit ), at which condensation in the 
second component ceases, varies with the discrepancy pa- 
rameter a = Nl ^^ 2 ■ The variation of (H) crit , nondimen- 
sionalised with respect to the value at a = ((H) { ), 
is shown in the inset of Fig. [TJ The figure clearly illus- 
trates that (J? C rit) deviates significantly from its value at 
a = as the discrepancy parameter is increased giving 
rise to a range of energies where only the first component 
condenses. 

In summary, we have derived a theoretical formulation 
of the thermodynamic state governed by the coupled NLS 
equations for a classical system of two weakly interacting 
waves in 3D. Numerical simulations are in excellent quan- 
titative agreement with the theory that is based on the 
equilibrium solutions of the kinetic equations for the NLS 
system of equations. The study presented here is relevant 
in quantifying condensation in a number of physical sys- 
tems. Notable examples include mixtures of BEC gases 
at finite temperature and nonlinear optics provided many 
modes are present that can be modeled semi-classically. 
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many useful discussions. 
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FIG. 1: Mass fractions as a function of averaged total energy 
for -^L = 0.5, -^r = 0.25, a = 0.1; inset - non-dimensionalised 
critical energy as a function of discrepancy parameter. 
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